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^-H Abstract 

o 

Source localization by matched-field processing (MFP) generally involves solving a 
^ ^ number of computationally intensive partial differential equations. This paper intro- 

^ duces a technique that mitigates this computational workload by "compressing" these 

X/^ computations. Drawing on key concepts from the recently developed field of com- 

I pressed sensing, it shows how a low-dimensional proxy for the Green's function can 

be constructed by backpropagating a small set of random receiver vectors. Then, the 
^ ^ source can be located by performing a number of "short" correlations between this 

P i proxy and the projection of the recorded acoustic data in the compressed space. Nu- 

^ merical experiments in a Pekeris ocean waveguide are presented which demonstrate that 

Q this compressed version of MFP is as effective as traditional MFP even when the com- 

pression is significant. The results are particularly promising in the broadband regime 
I where using as few as two random backpropagations per frequency performs almost as 

^ well as the traditional broadband MFP, but with the added benefit of generic appli- 

00 cability. That is, the computationally intensive backpropagations may be computed 

offline independently from the received signals, and may be reused to locate any source 
within the search grid area. 

ON 

2 1 Introduction 

J> 1.1 Background and Motivation 

X 

Matched field processing (MFP) continues to serve as one of the most widely used methods 
for localizing undersea targets acoustically. However, as the models governing undersea 
acoustic interactions become more sophisticated, often requiring fine-grain solutions to 
more complex partial differential equations, the tradeoff between run time and performance 
begins to worsen, perhaps unnecessarily. We will begin by discussing why this is the case 
and giving an overview of our approach to mitigate the problem. 

MFP generalizes standard array beamforming methods (e.g. plane wave beamforming) for 
locating an acoustic source in a complex environment (such as a multipath shallow water 
waveguide). MFP has been studied extensively both theoretically and experimentally as 
described in several review articles |25| EJ [3j . MFP is usually implemented by systematically 
placing a test point source at each point of a spatial search grid of L candidate locations, 
computing the acoustic field (replicas) at all the elements of the receiver array and then 
correlating this modeled field with the data from the real point source whose localization 
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is unknown to determine the best-fit location (see Fig. [T]). This approach works well when 
that the computational replica environment is sufficiently accurate. However, this direct 
implementation of MFP using brute force search would require L computation runs which 
can become numerically cumbersome for large search space especially when simulating 
complex propagation environments. 

One alternative to this direct implementation of MFP is to use a "backpropagation" algo- 
rithm (also referred to as "time-reversal imaging") to locate the unknown source. In this 
case, a time-reversed version of the recorded data is used as an initial waveform excitation 
along the array aperture using the principle of superposition, and then subsequently "back- 
propagated" numerically in the replica environment towards the grid search area [T7] . The 
unknown source location is then estimated from the maximum of the distribution of the 
backpropagated peak amplitude (or energy) across the grid search. Consequently, when 
compared to the direct implementation ffist mentioned, this backpropagation approach 
appears attractive at ffist glance, since it requires one computational run per unknown 
source. Nevertheless, this backpropagation approach becomes computationally expensive 
if multiple sources need to be located repetitively over the same search grid as the number 
of required computational runs would grow proportionally. For instance, this may occur 
when one tries to locate a source moving along a long track throughout the search space. 
Indeed, in order to be able locate any source throughout the search space using receivers, 
MFP would require computing N backpropagations by using sequentially each individual 
receiver as a backpropagation source |25| [2]. This would allow determining the full set 
of Green's functions associated with the channel between each search location and each 
receiver element. Alternatively, one could weight spatially the amplitude of the backprop- 
agated signals along the receiver array using N different orthogonal codes (e.g. obtained 
from an Hadamard basis). 

This article develops instead a compressive MFP formulation which reduces this compu- 
tational burden by pre-computing the backpropagation of a number M ^ of random 
test signals. The results of these backpropagations effectively encode the Green's function 
associated with the channel, and they can be re-used in subsequent localizations without 
any additional computational cost. This approach is inspired by recent work in the field 
of compressed sensing [QlTllIS], whose central message is that random projections provide 
an effective encoding for sparse signals. The motivation for compressed sensing is typically 
concerned with reducing the cost of acquiring signals by shifting the workload from sensor 
hardware to software ^i24j|20j, and is natural in applications where physical measurements 
are expensive compared to numerical computations. Here we explore a variation on this 
theme: mitigating the computational workload in software instead of the sensing workload 
in hardware. The proposed compressive MFP allows us to estimate the underlying ambigu- 
ity function central to conventional MFP algorithms over the entire search space using only 
M computational runs instead of A^, an effective speedup of a factor of N/M. In practice, 
these M simulations can be independently computed as a background process offline before 
the actual source signal is received. 

1.2 Related Work 

In this paper, we effectively demonstrate how classical localization procedures under a least- 
squares framework such as matched-field processing (MFP) may be solved in a reduced- 
dimensional space even without a-priori knowledge of the "best" dimension-reducing trans- 
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form. This property has been shown in similar forms in the mainstream canon of Com- 
pressed Sensing (CS) hterature. Davenport et al. have described a number of useful vari- 
ations on the theme of CS including a matched filtering detector. They have also 
described the "smashed filter" that is designed primarily for classification between a finite 
number of sets, but could easily be extended to parametric estimation. Wakin has also 
established some rigorous results on parameter estimation that relate the recovery prop- 
erties of a general compressive estimation problem to the properties of the manifold that 
these parameters induce. This work could be used to analyze this problem via its manifold 
parameters |26j . 

Carin et al. have utilized CS principles to show how a Green's function of a scattering 
field that is compressible in the wavelet domain may be recovered from a small set of 
measurements, though they use incoherence in their structured measurements to recover 
a scattering field, while we primarily care about the location of the source |T0]. Likewise, 
Marengo et al. have applied compressed measurements to the scattering problem, utilizing 
the target-sparse model to improve their performance |23j . 

Our work may also be viewed in the context of randomized SVDs [19! . ^^^^ field of 
research, the idea is to apply the matrix A to a series of random vectors $m as 
in order to determine the range space of A. For example, Chaillat et al. show how the 
inverse medium problem can be simplified using a dimension reducing random projection 
and solving the inverse problem in the reduced range-space [1^. Similarly to this field, we 
apply the time-reversal or adjoint of the Green's function Guj to random vectors in order 
to discover the range space of admissible ambiguity functions. 

There is also a large amount of recent research performing multi-target tracking under the 
"target-sparse" assumption. That is, the methods propose to simultaneously localize several 
targets that lie on some grid (or generally some set of points) by solving an ii minimization 
program. The recovered support resulting from this optimization corresponds to the grid 
points that the various targets are estimated to occupy. All of this work dovetails in very 
nicely with the main results of Compressed Sensing, which can be effectively leveraged 
to prove that the targets may be perfectly localized with high probability. Often, the 
painstaking effort in these papers involves showing that the Restricted Isometry Property 
(RIP) holds for the observation matrix. For example, Fannjiang et al. show the conditions 
under which a sufficiently small coherence is achieved for perfect recovery |16| . Gurbuz et 
al. show similar results for a Compressive beamformer, requiring a number of measurements 
on the order of the number of sources [18], but the application there is different in that 
they utilize a signal common to all sensors with an unknown time shift to localize their 
target in angle (assuming free space propagation), and apply the compression operator in 
time per-sensor instead of applying the operator across the range of sensors as we do. Also 
from a communications perspective, Cevher et al. demonstrate the relatively low amount 
of information to be transmitted for purposes of localization when using a Compressed 
Sensing framework [11] . These "target-sparse" approaches depend on targets lying exactly 
on the grid points. Also, by necessity these grid points must be spaced sufficiently far away 
from one another to avoid coherence-inducing correlations in the observation matrix. This 
creates a restrictive model of limited applicability. When a target is somewhere in between 
a set of grid points, the necessary conditions for recovery may not even approximately hold, 
similarly to how a discrete sinusoid corresponding to an off-grid point in the DFT will not 
be sparse in the frequency domain (or any other basis for any standard transform for that 
matter) due to DFT leakage. In contrast to this approach, we do not require our target 
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Figure 1: Schematic of a matched-field processing implementation in an ocean waveguide. The 
signal transmitted by a source (star symbol) located at an unknown location tq is recorded along 
a N elements receiver array after multipath propagation. Using a computational model of the 
original ocean waveguide, the location fo may be inferred by matching the actual received signals 
with the simulated replica waveforms obtained from varying the test source location (dot symbols) 
f throughout the search grid area. 

to lie on a grid point. However, instead of promising perfect recovery, we instead content 
ourselves to claim that our target may be localized to within a small neighborhood of the 
actual source location, or at least the location found via deterministic means. 



1.3 Outline 

The remainder of the paper is organized as follows. Section [2] briefly describes conventional 
MFP formulation for locating both single-frequency (narrowband) and broadband sources. 
Section presents the corresponding compressive MFP (cMFP) formulation for both cases. 
Section H presents numerical simulations in a Pekeris waveguide |21t pg 540-552] illustrat- 
ing the performance of cMFP in comparison to the conventional MFP results including the 
effects of additive ambient noise to the data and model mismatch due to uncertain knowl- 
edge of the actual environment. Section [S] extends this compressive approach to adaptive 
MFP. Section [6] summarizes the findings and conclusions drawn from this study. 



2 Conventional MFP 

A brief summary of the conventional MFP formulation is presented hereafter based on the 
standard solution of the linearized wave equation. The acoustic pressure field y{r, t) at a 
fixed point r and time t produced by a point source located at tq satisfies: 

1 dMf,t) _^2y^_,^^^^^^^^^^_._ 



^(f) 



where c(r) is the speed of sound and a(t) is the signal emitted by the source. The time- 
domain Green's function for the same environment g{'r,fo,t) is, by definition, the solution 
of Eq. ([l]) for a impulsive point source (i.e. for a{t) = S{t)) that satisfies all boundary 



conditions [211 Pg 540-552]. Using Eq. ([T]) (and assuming that the radiation condition 
apphes as ||r|| — )■ oo) the Fourier transform of the recorded pressure field at fvi, the n*'' 
element of a receiver array (n = 1..N) (see Fig. [T]), is denoted yuj{fn) and given by: 

yw{rn) = auiOujirn, fo) (2) 

where u is the frequency. The variables and gaj(^n>n)) denote respectively the Fourier 
transform of the source signal and time-domain Green's function. Using vector notation, 
Eq. (§ can be restated as: 

= a^G^{n)), (3) 

where is a (A'^ x 1) column vector obtained by stacking the complex amplitudes yaj(?n) 
measured along the receiver array. Similarly, the (A^ x 1) column vector Gcj{r) contains 
Green's functions guj{rn,^ between the receiver array elements and a source located at 
tq. Note that the position vectors are written in lowercase letters with arrows and the 
column vectors are written with capital letters in the remainder of this article. 



2.1 Single- Frequency MFP 

We start by considering the simplest MFP that works from measurements at a single 
frequency a; (as in ([s])), known as the harmonic (or narrowband) formulation. Given a 
set of measurements Y^j G across the N receivers at frequency w, we search for the 
location r in our region of interest TZ (and complex source amplitude /?) that best accounts 
for these measurements by solving the least-squares problem 

argminmin \\Yuj — BGuj(r)\\'^ ■ (4) 
ren /3ec " ^ ^ ^ 

With the location r fixed, the inner optimization problem is simply finding the closest point 
on the line spanned by G^^{r) to the point 1^. Plugging in the closed- form solution to this 
problem (see Appendix \K\ , the problem above reduces to: 

argmm ,, = argmax , (5) 

^ \\Guj{r)Y f \\G^(r)\\^ 

(where Y^ denotes the Hermitian transpose) which we will refer to as the normalized am- 
biguity function, and will refer to its maximization as normalized Matched Field Processing 
(nMFP). We show an example of the normalized ambiguity function in Fig. [2ja. 

The term Y^Guj{r) can be computed at every location r in an efficient manner using 
time- reversal. Precise values for ||Gtj(r)|p are typically not available when computing the 
backpropagation y^G^(r). However it is often the case (and we will assume this here) that 
these energies either do not vary much across our locations of interest, or vary predictably 
(e.g. cylindrical spreading of the field amplitude). Dropping the denominator yields the 
so-called unnormalized ambiguity function (alternatively the unnormalized Bartlett formu- 
lation) \25\ [2l [3], the objective function used for estimating the source location: 

r = argmax I /i(r)p where h{r) = Y^ Guj{r), (6) 
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2.2 Broadband MFP 



Now suppose that most of the energy of the source signal occupies some continuous band- 
width [wmin f^max], known as the broadband formulation. Ideally, we would solve Q over 
a continuum of uj values. However, for the sake of source localization, it is computationally 
advantageous to sample this bandwidth at K frequencies oji,oj2, - ■ ■ ,ojk, yielding K mea- 
surement vectors 1^^, where k G {1,2, ...K}. In this way, we can achieve a computational 
complexity at most K times the single frequency case, without sacrificing much precision. 

We now search for the location r that jointly matches the joint behavior of the measurements 
over multiple frequencies uji,i02, - ■ ■ ,i^K- The least-squares problem from Q becomes 



K 



arg mm mm 



(7) 



fc=i 



The inner optimization problem is separable over the (3^^, , and so the above is equivalent to 

K 



arg min V min 1 1 Y^^ - /3^^ G^^ (f ) 



k=l 



fc=l 



(8) 



As before, if the energies || 0(^^(0 IP homogenous across space and frequency, then a 
reasonable unnormalized approximation to the above is 



K 



K 



argmax^ I yj^G<^^(r) 1 2 = ^\h^^{i 



(9) 



k=l 



k=l 



The formulation in ^ assumes that the source amplitudes jS^^f, are unknown. If we have 
knowledge of the source signal's complex amplitudes, that is we know them up to a common 
amplitude and phase, then Q can be refined to 



arg mm mm 

r /3eC 





Yij^ 




aui-i_Guii (r) 






Y 


















Y 









(10) 



where the source amplitudes a^jj. are fixed and known. Applying again the results from 
Appendix |A| the inner optimization program can be solved in closed form, and so (10) is 
equivalent to 



arg max 



2^1.-1 la, 



k=i \^i^k\^\\Gu^dr)\\^ 



(11) 



as shown in Fig. [2]b, which we can approximate (by removing the denominator) as its 
unnormalized counterpart 

K 2 

'^^k ^^k 

(f) . (12) 

k=l 

Hereafter, we will refer to ([8| and Q as the incoherent MFP formulation, and (11) and 
(12) as the coherent MFP formulation. 



arg max 
f 
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Figure 2: Single- frequency (left column) and broadband-coherent (right column) ambiguity func- 
tions. These ambiguity functions shown on the dB scale ('201og]^Q(-)j for: (a, b) the standard MFP 
as described in Eqs. ^ and ( |TI| ), and (c, d) cMFP as described in Eqs. ( [T5| and (22 1 for the 
single-frequency case with M — 10 and broadband coherent case with M — 2 measurements per 
frequency, and (e, f) cMFP for the single-frequency case with M — 30 and broadband coherent case 
with M = 20 measurements per frequency. 
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3 Compressive MFP 



In this section, we describe Compressive Matched Field Processing (cMFP). This is an effi- 
cient method for acquiring a compressed version of the Green's function operator Gui{r) that 
exhibits a behavior in some regards similar to the dimension-reduced counterpart achieved 
via Principal Component Analysis, but may be obtained with only incomplete knowledge 
of the Green's function Gi^{r). Our approach works by precomputing the backpropagation 
of a small number of hypothetical received signals to construct a dimension-reduced proxy 
for the Green's function. Then, given the actual observed data 1^ we localize the source by 
finding the closest match between the received signal and the Green's function in the com- 
pressed domain. With the compressed version of Gi^{r) in hand, locating the source only 
requires computing a series of short inner products. In addition, the compressed version 
of Gu){r) is independent of the received signal, and so can be pre-computed and re- used 
for later observations. As we will demonstrate in Section |4j this cMFP strategy is effective 
even when the number of pre-computed random backpropagationsis far fewer than what 
would be required for a full acquisition of Gco{r) over the whole search grid area. 



3.1 Single- Frequency cMFP 



We start by discussing the single-frequency case in detail. First, we compute the compressed 
Green's function ^Guj{r), where ^ is a, M x N encoding matrix. Note that matrices 
are written in boldface letters in the remainder of this article. We construct ^Gu]{r) by 
backpropagating (i.e. applying G^^ to) a series of test vectors 'I'l, . . . , £ — we will 
discuss how the are chosen in the next section. 

The result of one of these computations ^!^Gi^{r) is a complex- valued acoustic field over f 
and requires as much effort to compute as the ambiguity function h{r). We stack up the 
results of these precomputations as rows in the ensemble 



[<i>i $2 



(13) 



This ensemble gives us access to an indirect, dimension-reduced version of Guj{r)- 

Given observations , we search for the r that best explains these random backpropaga- 
tionsin the compressed space. The least-squares program ^ becomes 

/3*G^(r)f, 



arg mm mm 

r (3 



(14) 



which, again using the results from Appendix [A] reduces to 



(narrowband cMFP) arg max 



where h{r) 



(15) 

The function h{r) is shown in Fig. ^c and ^e, and can be interpreted as a compressed 
version of the ambiguity function h{r) in ([g]) shown in Fig. [2ja. The cross sections in range 
and depth of these ambiguity functions are shown in Fig. [3|a and[3]b. 
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Note that unlike the standard MFP, in this case the precomputations give us direct access 
to the denominator ||$Gtj(r)|p (we simply take the norms of the columns of $G(^(r)), 
and so we leave it in the optimization program. As shown in the results section, this 
normalization term plays an important role in improving the source location estimation 
when the magnitude of the Green's function varies significantly across the search grid area. 



Notice that an evaluation of (15) at a point r essentially only requires an inner product 
between the encoded observations and the M- vector $G(^(r) formed from the back- 
propagated fields at point r from all M test vectors. 



3.2 Random Projections 

The question remains as to how to choose the encoding matrix $ so that solution to 
the cMFP (15) is the same (or close to) the solution to the standard MFP The 



corresponding least-squares problems are 

standard MFP : argmin^^/j \\Y^ - /3G<^(f)||2 (16) 

cMFP : argmin,-,;3 ||* (F^ - /3G^(r)) f • (17) 

These two programs will have similar solutions if their functionals are close to one another 
for all values of /3 and r. If Y^j = aGt^(ro), then the performance of the cMFP will match 
that of the standard MFP when $ preserves the energy of the differences between the 
observations 1^ and all scalar multiples of the Green's function at different points: 

\\^{Fi - F2)f ^ \\Fi - F2f for all Fi, F2 e F := {F : F = aGUr), a e C; r £ U}. 

(18) 

Essentially, we want $ to stably embed (i.e. preserve the distances between members of) 
the set F into C^^. 

We propose taking $ to be a random linear mapping. This choice is inspired both by 
classical results in theoretical computer science and from the recently developed theory of 
compressive sensing. In the mid-1980s, Johnson and Lindenstrauss [22] demonstrated that 
the distances within a finite set of n points are essentially preserved through a random 
projection into a space of dimension ~ logn (see also |13l [T]). Recently it has been shown 
that this same type of projection also embeds sparse signals into a low-dimensional subspace 
[3], a result which plays a key role in compressive sampling [TJ [9], and are effective at 
reducing the dimensionality of certain types of manifolds [S]. 

We will discuss the particular the case where $ is a random orthoprojection, although the 
results will be almost identical for many different choices of random $ (e.g. with entries 
that are independent and identically distributed Gaussian or ±1 random variables). To 
generate we simply draw an Af x matrix of independent Gaussian random variables 
with unit variance, orthonormalize the rows using the Gram-Schmidt (or QR) algorithm, 
and then multiply by N/M. For an arbitrary fixed vector F, the random orthoprojection 
$ obeys two properties ^3j: 

E = \\Ff (19) 

P {| - \\Ff I > e} < 2e"^ (20) 

This allows us to interpret the compressed energy functional ||^(5^ — f3G^{r))\\'^ in ( [l7[ ) as 
a random process, indexed by (3 and r, whose mean is the standard energy functional ||K<j — 
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Figure 3: Cross-sections of the ambiguity functions displayed on Fig. (a, b) single-frequency 
case described by Eqs. ^ and Eqs. (151; (c, d) broadband coherent case Eqs. (11) and (22); range 
(left column) and depth (right column). Here we show the normalized standard MFP (nMFP) and 
the cMFP (cMFP) for various values of M. The dashed lines show the boundaries for the main lobe 
and region of uncertainty that we are able to localize within under the presence of modest noise. 
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PGu]{r)\\'^ in ( 16). At a fixed point j3, r, this random process is concentrated around its mean 



roughly hke a Gaussian random variable with standard deviation ^^2/M\\Y^^ — /3Ga;(r)||. 
The larger we make M (the more random vectors we precompute backpropagations for), the 
tighter the concentration. By construction, when M = N, = I and we have acquired 

a "lossless" version of the Green's function G^{r)^ meaning that the functionals are exactly 
equal to one another. In general, however, we will be interested in cases where there is 
a significant compression factor M <^ N and benefit from the associated computational 
savings. 



3.3 Broadband cMFP 

The cMFP formulation can be readily extended to combine observations at multiple frequen- 
cies in both the incoherent and coherent cases. For frequencies wi, W2, . . • , ^j^k, we generate 
a sequence of M x random matrices , ^ui2 1 • • • ) ^ojk backpropagate the rows 
of each (for a total of MK time-reversals) to acquire ^ujiGoji{r), . . . ^^ujkGujk{^- Then 
given observations Y^-^ , . . . , Y^jj^ , we compress them by calculating , ■ ■ ■ , ^uik ' 

and then using the compressed versions of the G^^^. , we proceed as in ([7]) for the incoherent 

(incoherent cMFP) ar, 

= arg min ^vca^W Y^^ - /3^, G^^ (f ) 



case 



K 



gmin min ^W^^j^Y^k - P^k^^kGuk{r)\? 



K 



k=l 
K 



^,yj^*^^*.,G.,(f)|2 

arg max \ 



k=l 
K 



arg max 



E 



^ \\^u}kGuj,.{r) 



and as in (10) for the coherent case: 



(21) 



(coherent cMFP) arg min min 











s> Y 






s> Y 





(^UJl ^(jJl GuJi ij") 
^U)2 G(jj2 (^) 



P^UJk ^'^K GulK 



arg max 



Efc=i otu, X., G^fe (f ) 



E 



K 

k=l 



arg max 



E 



K 
k=l 



a, 



^^kG^kim' 



(22) 



The incoherent and coherent respectively illustrated in Fig. [2jd and [2jf and in 

Fig.[3jc and[3jd. Note that in this coherent case, the optimization is identical in its structure 
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to the single-frequency case. In particular, by concatenating: 









'Y 






G(f) = 


0/.UJ2 GuJ2 (^) 


Y = 


Y 


* = 


^L02 








Y 







(23) 



we see that the coherent broadband formulation (22) shares the same formulation as the 
single frequency case (15). 



4 Numerical simulations 

In this section, we present numerical experiments demonstrating that underwater acoustic 
sources can be localized from highly compressed versions of the Green's functions. Our 
cMFP results give locations estimates for single-frequency, incoherent broadband, and co- 
herent broadband that are comparable with the traditional MFP. After the initial pre- 
computation (which consists of backpropagating the random codes at each frequency), the 
cMFP is substantially faster than the traditional MFP, requiring only a short inner product 
to be calculated at each search location. 

The MATLAB code generating all the numerical results presented in this section is available 
online 



4.1 Numerical set-up 

All numerical simulations were conducted using a 200m deep Pekeris waveguide and the 
Green's functions were computed using a standard normal mode code [21\ pg 540-552]. 
The two dimensional search grid area in depth and range spans respectively [10m 190m], 
and [5000m 5810m] for the single frequency and broadband incoherent simulations. The 
range span for the broadband coherent simulations was reduced to [5000m 5270m] to keep 
constant the number of search locations over which the ambiguity functions are computed 
since the effective resolution of the ambiguity function in the coherent case was about 3 
times higher in range (see Fig. [2] and Fig. [3]). A uniformly spaced vertical line array with 
= 37 elements spaced between 10 and 190 meters was used to sample the acoustic field. 
The Green's functions between each of the search locations and the receiver array (see Fig. 
[T]) were calculated across K = 20 different frequencies between 141 Hz and 160 Hz (the 
narrowband configuration uses 150 Hz). Given the selected numerical set-up, the natural 
resolution in frequency of the computed Green's function is around 5 Hz; that is, Gtji(r) 
and Gi^2{r) are essentially uncorrelated when |a;i —uj2\ > IOtt. The selected sample spacing 
of 1 Hz falls well within this frequency resolution. 

After selected a source location ro inside the region of interest, observations at the K 
frequencies were simulated using the forward model, and uncorrelated zero-mean Gaussian 
noise was added to the result: 

Y^, = a^^G^^ifo) + Zk, ZfeGC^, Zfc~Normal(0,a2/), (24) 



^Download the code at littp://users. ece.gatech.edu/~wmantzel3/cmfp/code. zip 
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Figure 4: DeEnition of the elliptical distance metric used for the performance study of cMFP. 
Because range errors tend to be greater than depth errors in long-range localization estimates, 
our results use an elongated distance metric that gives greater weight to the depth than the range, 
leading to unit balls that are ellipses instead of circles as shown here (see Eq. ( p7| ). The color scheme 
for this ambiguity surface 101og;^Q(|/i(r)p) = 101ogjg(|l^G'^(r)p) has been lightened somewhat to 
allow for better visibility for the overlaid ellipses. 



where each Zk has i.i.d. Gaussian real and imaginary parts with variance In all of 

our experiments, we set a^j^ = 1 for all k. The signal-to-noise ratio (SNR) corresponding 



to noise variance a"^ is 



SNR= 101og.„( l°"l'l;^;f°>"' ) (25) 

in the single frequency case, and 

SNR = lOlogio |^ ^^=i'"-^'^|^^'-'-^^°)"' ^ (26) 
in the broadband case. Unless otherwise stated, we used an SNR of 16 dB. 



Given a set of observations, we estimate the source location by solving (15) single fre- 
quency), (21) (broadband incoherent), or (22) (broadband coherent) and compare against 
the standard MFP formulations ([6]), and (12) As stated, these optimizations problems 
are over a continuous variable r — in practice, we compute these functionals on a finite 
grid of points and choose the maximum from amongst these points. We used a 90 x 90 
grid for the simulations presented below, which corresponds to 2m spacing in depth, a 9m 
spacing in range in the single-frequency and broadband incoherent cases, and 3m spacing 
in range in the broadband coherent case. We wish to emphasize that while our solution 
will of course lie on one of these grid points, the actual source location is chosen to be an 
arbitrary point. 

The natural resolutions in depth and range of the ambiguity function h^{r) differ, as shown 
in Fig. |4] and Fig. [s] for a source located at ttq = (5540m, 100m) in (range, depth) for 
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a single frequency u = SOOvr rad/sec (150 Hz) for a source located at = (5540, 100) 
in (range, depth). In this case, the main lobe has a width of ~ 360 m in range and ~ 
32 m in depth. Again the grid spacing of 9m /2m in range/depth falls well within this 
resolution. The spatial resolution of the ambiguity surface in the selected multi-modal 
Pekeris waveguide is primarily a function of the source-receiver array configuration as well 
as the selected frequency band [25, 2j In light of these differing spatial resolutions, we use 
a weighted norm to report distance errors in most cases presented here in this section. The 
distance from a point tq = (rg^"^*^, Tq^^**^) to the estimated source location r is computed 
using the elliptical distance: 



range -range \ 2 /depth -depth \ 

W r° .Hel (27) 



^ grange J ^ gdepth 



We use e^'^P*'^ = 3m and e'''^^^^ = 36m for the single-frequency and incoherent cases, and 
grange _ 12^^ [j^ coherent case. The values of e'^^^^^ and e'"'^"^'^ were chosen so that 
the contour {r: Ht'o — r||e = 1} was approximately the same as the isosurface of the 
ambiguity function at 0.9 of its maximum. Equidistant points from ro = (5540, 100) for 
11?^ — r||e = 1, 5, and 10 are shown in Fig. |4j For example, an error of 14.4 meters in range 
and 0.9 meters in depth translates to 0.5 units of distance error in the elleptical || • ||e norm. 



4.2 Localization performance of cMFP. 

Fig. [5^ compares the performance of cMFP (see Eq. ( 15 )) and MFP (see Eq. (eq:amb-norm- 
eq:amb)) for locating a harmonic source (/ = 150Hz). The SNR of the received data vector 



(see Eq. 24) was set to 16 dB. For a fixed M we aggregate performance statistics across 
1000 simulations: 100 different source locations (chosen from TZ uniformly at random) and 
10 different draws of the for each location. For each test simulation, the error between 
the true and estimated target location was recorded in units of the target ellipse radius (see 
Fig. [4]). From the results of the 1000 simulations, we calculated the empirical distance tail 
probability P/vf ((i)-for a given number of random backpropagations M- as the fraction of 
results that produces a location estimate r with |[r — ro||c > d. As shown, we are able to 
estimate the target within the unit ellipse more than 99% of the time from only M = 6 test 
vectors. Notice that the cMFP actually outperforms the unnormalized version of the MFP 
(from ^ above) when M ~ 6. This happens because the cMFP has an estimate of the 



normalizing factor in the denominator, as shown in (15). The cMFP is really an estimate 
of the normalized MFP in ([s]), and indeed that formulation is what the cMFP approaches 
as the number of random backpropagations M becomes equal to number of receivers A^. 

The cMFP was also tested in a variety of SNR for the single- frequency case. Fig. [Sja shows 
the probability that the localization estimate is within the first ellipse (i.e. d < 1) as a 
function of the number of random backpropagations M. In all cases, the failure proba- 
bility asymptotically decreases exponentially in the number of random backpropagations. 
Finally, Fig. ^ shows the tail probability of distance error for a fixed number of random 
backpropagations M = 20. As expected, the performance of cMFP gradually decrease as 
the SNR of the measurements is reduced from 16dB to OdB, similarly to what occurs when 
using conventional MFP |3]. 

Fig.[6]and Fig.[7]show a similar performance study for respectively the broadband incoherent 
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Figure 5: (a) Tail probability of distance error \\r — ro||e (see Eq. (27)) for the single-frequency 



cMFP formulation (see Eq. (15)) at 150 Hz. Pkiid) is the probability that the localization is worse 
than some distance d using M random backpropagations. The dashed lines indicate the performance 
under normalized and unnormalized MFP (Eq. ([5| and Eq. Qj. The next two plots show results 
for Pnid) over various SNRs of the received signal with (b) fixing d — 1 and (c) fixing M — 20. 
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Figure 6: Same as Fig. [s] but using instead the incoherent broadband cMFP formulation (see Eq. 
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Figure 7: Same as Fig. [s] but using instead the coherent broadband cMFP formulation (see Eq. 
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cMFP (see Eq. (21)) or broadband coherent cMFP(see Eq. (22)) formulations, including 
the influence of the SNR of the measurements as well as the number the number of random 
backpropagations M. Note that in Fig. [7| the horizontal axis is normalized differently than 
in the other two cases due to the different spatial resolutions of the ambiguity surfaces 



(see Eq. (27)). Our intentions are not to directly compare the performance of each of the 
three cMFP formulations (the coherent localization being always better as expected), but 
rather to show that in each case the selected cMFP formulation performs as well as the 
corresponding normalized MFP formulation and better the corresponding unnormalized 
MFP formulation. This is especially true for the broadband coherent cMFP results as 
Fig. [7ja shows that with just M = 1 measurement per frequency, we achieve an error 
within 3 times what standard MFP gives us at least 90% of the time, and with M = 2, we 
fall within about 10% distance error of what MFP gives us about 99% of the time. 

Furthermore, note that we do not show results for M = 1 for the broadband incoherent 



cMFP formulation (Fig. [oja) as in this case ^kGuj^ (r) is a scalar for each r , and ( 21 ) reduces 
to 

arg max > tt^—^ — t^ttt^ = arg max > r^— — -—p, (28) 

r ||$fcG^,(r)||2 ^ r ^ |*fcG^,(r)|2 

K 

= argmaxV|$fcK;j2_ ^29) 

r ^ — ' 

k=l 

This optimization problem is ill-defined, as the functional does not depend on r. 



4.3 Evolution of the main lobe to side lobe ratio of the cMFP ambiguity 
surface. 

Fig. [8] shows the logarithmic variations of the main lobe to side lobe ratio of the ambiguity 
surface obtained with the single frequency and broadband coherent cMFP formulations for 
increasing number of random backpropagations M. In each case, the displayed values rep- 
resent the median value of the main lobe to side lobe ratios obtained from 1000 simulations 
for each value of M. Here the main lobe is defined as the maximum of the ambiguity surface 
\h{r)\ (obtained from the corresponding conventional MFP formulation, e.g. see Fig. [2^-b 
and Fig. [3]) over the region of interest IZ, and the side lobe as the maximum of \h{r)\ over 
the search area TZ excluding an ellipse E of the approximate size of the main lobe. We 
show the cross sections of the ambiguity function in Fig. [3] where we illustrate our choice of 
main lobe ellipse parameters that define our main lobe ellipse E. For the single frequency 
case, the ellipse has parameters e'''^"^'^ = 180 meters and e'^^P*'^ = 16 meters (the broadband 
coherent case uses e^'"^^^^ = 72 meters and e^'^P*'^ = 16 meters) as illustrated in Fig. [sj The 
logarithmic value of the main lobe to side lobe ratio is computed as: 

201og.„(^=?i^^), (30) 

Vmax^g7^\£; \h{r)\J 

Note that for small M values, the cMFP side lobes may be significantly larger than their 
standard MFP counterparts. The concentration inequality ( |20[ ) suggests that as M gets 
larger, the side lobes dampen. This behavior is observed in Fig. [8) Note that since the $ 
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Figure 8: Evolution of the main lobe to side lobe ratio (in dB) of the estimated ambiguity surface 
(e.g. see Fig. ^ vs. number of random backpropagations M using either (a) the single frequency 
cMFP formulation at 150 Hz or (b) the broadband coherent MFP formulation (see Eq. (22)). Note 
that in each case the main lobe to side lobe ratio of the ambiguity surface obtained with cMFP 
reaches the main lobe to side lobe ratio value obtained using the corresponding nMFP formulation 
(dashed line) when M = N = 37. 



matrix is an isometry when M = N, the side lobes in this case are exactly the same as for 
the standard MFR 



4.4 Influence of model mismatch on the cMFP performance 

Previous studies have shown extensively that one major liability of MFP is sensitivity to 
model mismatch which occurs when one has an incorrect model for the ocean waveguide (e.g. 
sound speed profile error) . Since MFP exploits the knowledge of the environment (via the 
Green's functions), its numerical accuracy must be sufficiently accurate, to ensure accurate 
source localization. Here we simply ensure that the localization accuracy of cMFP remains 
comparable to conventional MFP in the presence of error in the sound speed value. To do 
so, a set of received signals with a set SNR of 16 dB were computed for a reference sound 
speed of 1520 m/s. The broadband coherent cMFP -using M = 4 random backpropagations 



per frequency (see (22)) and normalized MFP formulation (see (|11[)) were then implemented 
using backpropagations in a simulated environment with different nominal values for the 
sounds speed (between 1520 m/s and 1530 m/s) than the reference value of 1520 m/s. Fig. [9] 
shows that the cMFP performs substantially the same as traditional MFP, for better or for 
worse. We show the average distance error in actual Euclidean distance (meters) as 



instead of ellipse distance. The small localization error occurring even without modeling 
error is due to the fact that the true source location did not coincide exactly with one the 
grid search location f. 

Note also that the range error tends to dominate for sufficiently large modeling error: the 
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Figure 9: Evolution of the localization error for broadband coherent cMFP and corresponding 
conventional MFP a for increasing error of the modeled sound speed value. The correct sound speed 
value is 1520m/ s here. Notice that the localization errors obtained from cMFP (circle symbols) 
match closely the localization errors obtained obtained from standard MFP (cross symbols). 



slope of the displayed error values is roughly 5000/1520 (the nominal range divided by the 
nominal speed of sound) as we would expect because a 15 m/s error in the speed of sound 
of 1520 m/s causes a corresponding approximate 1% distortion in the apparent range, or 
50 meters, i.e. 15/1520 • 5000. 



4.5 Application of cMFP for tracking a moving source. 

The advantage of cMFP over conventional MFP for locating a moving source along a long 
track is illustrated here. Fig. [To] displays the arbitrary path of a source moving along a 
parabolic trajectory (dashed lines). For the sake of simplicity, the Doppler effect is not 
accounted: this moving source scenario is simply simulated as 100 successive stationary 
sources located along the parabolic trajectory. For each positions, the SNR of the re- 



ceived signals at the vertical line array is constant and equal to 16 dB (Fig. 10 a) or 8 dB 
(Fig. |10[b). Conventional broadband coherent MFP is implemented by running 100 succes- 
sive backpropagations per frequency over the search grid to estimate the source trajectory 
(see crosshair symbols). On the other hand, broadband coherent cMFP is implemented 
using M = 2 random backpropagations per frequency to estimate the same source tra- 
jectory (see cross symbols). The median value of the distance errors (computed from Eq. 
(|31|)) between the estimated and actual source trajectory is Im when using both MFP and 



cMFP for a SNR of 16dB. A slightly higher error of 1.6m (resp. 1.1m) for the cMFP (resp. 



MFP) was found for a SNR of 8dB. Overall, Fig. 10 indicates that cMFP can potentially 
achieve comparable source tracking performance with a significantly reduced number of 
simulations. 
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Figure 10: Tracking of a source moving along a parabolic source trajectory (dashed line) using either 
coherent broadband cMFP, implemented with M = 2 random backpropagations per frequency for 
the whole search grid, or using conventional broadband coherent MFP. For each of the 100 source 
positions, the SNR of the received signals at the vertical line array is constant and equal to (a) 16dB 
or (b)8 dB. 

5 Extension to adaptive MFP 



Several variants of the MFP algorithm have been proposed in the existing hterature |25| 121] 
to enhance the robustness and performance of the basic Bartlett formulation presented 
above (see Eq. Q). This can be especially beneficial in the presence of added coherent 
noise to the received data vector 1^ (see Eq. ([s])). To do so, these higher resolution MFP 
algorithms are data adaptive, but typically have also a high resolution in their environmen- 
tal knowledge requirements. A commonly used adaptive MFP formulation is the Minimum 
Variance Distortionless Response (MVDR) formulation. The MVDR formulation adap- 
tively constructs a replica (or weighting) vector to yield a minimum mean square response 
to the recorded noise field along the receiver array while maintaining a constraint of unity 
processing gain for the incoming signal vector 1^ [21, pg 540-552]: 



\h. 



MVDR 



(f)p 



-1 



(32) 



where K is the N x N is the empirical correlation matrix from multiple realizations of the 
noisy received data vector Y^: 



(33) 



1=1 



The physical interpretation and performance analysis of the MVDR formulation (see Eq. 
([32|)) over the simple Bartlett formulation (see Eq. Q) have been discussed extensively in 
the previous literature |25l |2T] and thus will not be further repeated in this article. 

The previous cMFP formulation can be readily extended to handled adaptive variants of 
the simple Bartlett MFP algorithm as discussed in Section III.C. For instance, using Eq. 



(32) and by direct analogy to Eq. (15), the magnitude square of the compressive MVDR 



21 



ambiguity surface is: 



|/^WD«(^-.)|2 ^ (($G.(f))^ (*K*^)-' *G.(r))" . (34) 
So once we have computed the M tests measurements $Ga;(r), they can be readily apphed 



to either the compressive adaptive MFP formulation (see Eq. (34)) or the simple Bartlett 
formulation (see Eq. (flSj)) to locate the unknown source. 



6 Conclusions 



We have shown here how dimension-reducing random projections can greatly reduce the 
computational cost involved with source localization via matched-field processing. When 
compared to the location of the maximum of the ambiguity surface obtained from conven- 
tional MFP using N distributed receivers, the localization error achieved by cMFP scales 
down as square root of the number of random backpropagations M. The proposed cMFP 
formulation has also the added benefit to be able locate any source within the search grid 
area using only M random backpropagations, while conventional MFP would require at 
least backpropagations to do the same. Thus cMFP provides an effective speedup factor 
of N/M per frequency, which can be significant when a large number of receivers is avail- 
able to locate a broadband source. Consequently this cMFP technique enables the ability 
to both broaden the search space and employ more sophisticated models of the Green's 
function, without introducing worries about sacrificing real-time performance 

This compressive approach is not limited to source localization, and could be extended to 
a more general type of machine learning problem when matches are evaluated via inner 
products (or equivalently via Euclidean norms) . This type of approach has the potential to 
substantially decrease computational complexity in these cases, while admitting a negligibly 
small probability of error. 



A Closest point on a line 



For fixed vectors U,V £ C" , the following optimization program finds the closest point on 
the line spanned by f to 



min \\U-pvr. 

/3ec 



The functional above attains its minimum value of 

\\U 



2 l^^f/P 



when 

||y||2 

This fact can be verified by differentiating ||f7 — with respect to the real and imaginary 
parts of /3, and solving for value of f3 that makes them both equal to zero. 
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